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HIGHLIGHTS 


•  A  noble  multi-D  mathematical  model  of  large-scale  lithium-ion  batteries  is  developed. 

•  The  model  predicts  current,  voltage,  temperature,  SOC,  and  SOH  distribution. 

•  The  model  performs  simulation  under  dynamic  operating  conditions. 

•  Various  experimental  validations  confirm  that  the  model  has  high  accuracy. 
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Multi-dimensional  model  of  large-scale  lithium-ion  batteries  is  developed.  The  model  is  based  on 
equivalent  circuit  model  (ECM)  which  is  capable  of  dynamic  response  simulation.  Model  parameters  are 
functions  of  both  state  of  charge  and  temperature,  which  are  implemented  by  bilinear  interpolation 
method.  Local  degradation  effects  such  as  capacity  fade  and  resistance  growth  are  incorporated  into  the 
multi-D  model  by  empirical  method.  Local  distributions  of  current,  potential,  temperature,  state  of 
charge  (SOC),  and  state  of  health  (SOH)  are  observed  with  the  developed  model.  It  is  found  that  local 
degradation  worsens  under  severe  operating  condition.  Low  thermal  conductivity  of  a  battery  results  in 
non-uniform  SOH  distribution.  Simulation  result  shows  good  agreement  with  experiment  data  under 
various  operation  modes. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Since  General  Motors  successfully  launched  the  first  plug-in 
hybrid  electric  vehicle  (PHEV)  powered  by  LG  Chem’s  lithium-ion 
battery,  lithium-ion  battery  has  been  getting  more  public  interest 
as  a  strong  candidate  for  future  automotive  power  source. 
Compared  to  small  batteries  used  for  portable  electronics  which 
usually  operate  under  steady,  mild  operating  conditions,  automo¬ 
tive  batteries  are  exposed  to  more  severe  operating  conditions.  In 
addition,  automotive  batteries  must  sustain  at  least  ten  years  with 
guarantee.  For  these  reasons,  it  requires  much  time  and  efforts  to 
develop  lithium-ion  batteries  for  automotive  applications. 

Mathematical  modeling  approach  can  be  applied  to  battery 
engineering  to  save  time  and  cost  required  for  battery  develop¬ 
ment.  By  solving  governing  equations,  physics-based  model  [1-3] 
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can  predict  battery  performance  and  analyze  physical  behavior 
before  the  actual  battery  is  constructed.  When  a  battery  exists  and 
experimental  data  of  the  battery  is  available,  empirical  models 
which  mimic  the  battery  performance  can  be  generated.  Equivalent 
circuit  models  [4-7]  and  fitting  function  models  [8,9]  belong  to  this 
category.  As  empirical  models  are  relatively  simple  and  fast,  they 
are  usually  used  as  control  algorithms  for  battery  management 
system. 

When  the  size  of  the  battery  becomes  large,  it  is  important  to 
consider  dimensional  effects  such  as  battery  shape  and  lead-tab,  or 
the  high  voltage  outlet  position  on  local  distribution  of  current 
density,  electrical  potential,  state  of  charge,  and  state  of  health. 
Most  of  mathematical  models  listed  above  are  one-dimensional  or 
lumped-scale  models  which  are  inadequate  for  analyzing  dimen¬ 
sional  affects.  Therefore,  it  is  necessary  to  apply  multi-D  modeling 
to  design  large-scale  batteries.  In  development  of  an  automotive 
battery  pack  which  is  usually  cooled  by  liquid  coolant  or  air-flow, 
computational  fluid  dynamics  (CFD)  is  used  to  check  thermal- 
fluidic  behavior  of  the  battery  pack  under  various  operating  con¬ 
ditions.  In  order  to  accurately  predict  thermal-fluidic  behavior  of 
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Fig.  1.  Schematic  of  ECM  with  two  RC-pairs  representing  a  lithium-ion  cell  used  for  the 
present  study. 


the  battery  pack  under  given  load  condition,  it  is  important  to  know 
the  amount  of  heat  generation  which  is  directly  affected  by  dy¬ 
namic  battery  voltage.  This  battery  voltage  can  be  obtained  by  us¬ 
ing  multi-D  battery  model  incorporated  with  CFD.  Kwon  et  al.  [8] 
extended  Gu  [9]’s  circuit  network  based  Zn/NiOOH  cell  model  to 
two-dimensional  lithium-ion  cell  model  by  solving  positive  and 
negative  potential  field  instead  of  circuit  network.  This  kind  of 
model  is  typically  called  ‘NTGK  (Newman-Tiedemann-Gu-Kim) 
model’.  NTGK  model  solves  two  electric  potential  fields  coupled  by 
fitting  parameters  acquired  from  constant  discharging  profile  of  a 
cell.  This  model  is  useful  for  predicting  local  distribution  of  elec¬ 
trical  potential,  temperature  and  SOC  which  come  as  dimensional 
effects  in  a  large  size  battery.  However,  as  the  fitting  parameters  are 
functions  of  DOD  (depth  of  discharge)  and  they  do  not  have  any 
time-dependent,  relaxation  terms.  For  this  reason,  NTGK  model  has 
limitation  in  simulating  dynamic  response  to  transiently  changing 
load  profile. 

Mathematical  battery  models  can  be  applied  to  battery  life 
prediction  also.  Life  prediction  can  be  dealt  with  a  physics-based 
approach  which  requires  a  deep  understanding  of  degradation 
mechanism  taking  place  in  a  lithium  battery.  Actually,  this  degra¬ 
dation  mechanism  is  quite  complex  to  explain  with  physics-based 
governing  equation  yet.  Instead,  a  long-term  battery  life  is  typi¬ 
cally  estimated  with  empirical  methods  [10,11  .  As  a  large-scale 
battery  may  have  non-uniform  degradation  according  to  battery 
geometry,  material  properties,  and  operating  condition,  multi-D 
battery  model  can  be  useful  to  predict  local  degradation. 

In  this  paper,  we  introduce  how  to  develop  an  efficient  multi-D 
battery  model  which  can  be  applied  to  various  purposes  in  the 
development  of  large-scale  batteries. 

2.  Models 

2.1.  Equivalent  circuit  model 

Equivalent  circuit  model  (ECM)  expresses  a  cell  with  simplified 
electric  circuit.  Fig.  1  shows  a  schematic  of  ECM  with  one  serial 
resistance  and  two  RC  pairs  (five  ECM  parameters  total),  which  is 
used  as  a  baseline  mode  for  the  present  study.  Cell  voltage  can  be 
calculated  by  subtracting  voltage  losses  from  open  circuit  voltage 
(If)  which  is  a  function  of  SOC  as  shown  in  Fig.  2. 

Veen  =  u  -Vs  -V, 

r  u  =  /(soc) 

-V2  where  soc  =  soco__j_.  f  J(t).dt  (D 

l  time 


Fig.  2.  OCV  profile  of  the  lithium-ion  cell  used  for  the  present  study. 


From  Kirchhoffs  law,  three  voltage  loss  terms  (Vs,Viy2)  can  be 
found  as  follows: 


£  +  C2^  =  / 


Ri 


(2) 


Note  that  Vs  is  a  function  of  only  resistance  (Rs)  whereas  V\  and 
V2  are  functions  of  both  resistance  and  capacitance  which  has 
voltage  relaxation  effect  during  dynamic  operation  of  a  battery. 
Electric  current,  /  is  a  function  of  time,  not  a  constant  here. 

It  is  assumed  that  the  serial  resistance,  Rs  is  a  function  of  tem¬ 
perature  while  the  RC-pair  parameters  (Ri,R2,Ci,C2) are  functions  of 
both  temperature  and  SOC.  In  order  to  minimize  bias  from  tem¬ 
perature  rise  during  experiment,  ECM  parameters  are  acquired 
from  current-pulse  and  relaxation  test  (see  Fig.  3)  in  this  study. 
Iteratively  comparing  simulation  result  against  experiment  data, 
the  model  parameter  for  a  given  operating  condition  is  found  as  an 
optimal  fit.  However,  input  variables  such  as  SOC  and  temperature 
used  for  a  model  parameter  continuously  change  when  a  cell 
operates.  Therefore,  model  parameters  should  be  updated  at  each 
time  step  by  interpolation.  In  order  to  fit  two-dimensional  pa¬ 
rameters  such  as  RC-pairs  that  are  functions  of  both  temperature 
and  SOC,  a  bilinear  interpolation  method  is  used  as  follows: 

Table  1  is  a  two-dimensional  experimental  data  of  interested 
parameter  such  as  y  =  Ri(T,soc).  This  experiment  table  must  be 
filled  out  before  starting  interpolation.  Four  available  temperature 
(-7  °C,  0  °C,  25  °C,  45  °C)  cases  are  used  to  build  the  model  pa¬ 
rameters  for  the  present  study.  In  order  to  find  out  value  y  at  point 
(T,  soc)  where  experiment  data  does  not  exit,  first,  read  value  y  of 
four  neighboring  points. 


yu  =  R\(h,  soci) 
y\2  =  Rl(7l,  soc2) 
y2 i  =  R|CT2.  soct) 

y22  =  ^i(^2i  soc2) 

Then,  linear  variable  p  and  q  are  calculated  as  follows: 


P  = 


T-T! 


SOC  -  SOCi 
SOC2  -  SOC! 


(3) 


(4) 


t2-t i’ 
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Fig.  3.  Extraction  of  ECM  parameters  from  iterative  comparison  of  simulation  and 
experiment  data  (discharge  and  rest  test  @25  °C). 


Finally,  y  =  fti(T, soc)  is  found  out  by  plugging  Eqs.  (3)  and  (4) 
into  Eq.  (5). 


y  =  fii(T,soc) 

=  (i  -p){i-q)-yn  +P-0  -<j)-yi2  +  0  -p)qy2i  +pqy22 

(5) 

By  applying  the  above  bilinear  interpolation  method,  two- 
dimensional  ECM  parameters  extracted  from  the  experimental 
data  of  the  actual  Li-ion  cell  for  the  present  study  is  shown  in  Fig.  4. 
Various  different  cell  models  can  be  also  constructed  with  the  same 
manner. 


2.2.  Multi- dimensional  ECM 


Fig.  4.  ECM  parameters  acquired  from  bilinear  interpolation  for  the  present  study:  (a) 
R i  =J[ SOC,  T),  (b)  Ci  =J[ SOC,  T ). 


A  schematic  of  calculation  domains  in  the  present  multi-D 
lithium-ion  cell  model  is  presented  in  Fig.  5.  The  lithium-ion  cell 
model  has  three  domains  (cell  body,  positive  tab,  negative  tab). 
ECM  presented  above  is  inserted  into  each  mesh  of  the  cell  body 
domain  after  scaling  down  according  to  mesh  size.  In  order  to 
incorporate  all  of  small  ECMs  together,  the  following  electric  po¬ 
tential  equations  are  solved  in  the  cell  body  domain. 

'flit”  =  1  , } where  J  =  Jrxn  -Js  (6) 

<rnV z<Pn  =  -J  J 

In  Eq.  (6)  Jrxn  is  volumetric  electrochemical  reaction  current 
density  whereas  Js  is  volumetric  side-reaction  current  density 
defined  as  follows: 


Table  1 

Two-dimensional  data  table  for  bilinear  interpolation  of  ECM  parameters  for  the 
present  model. 


T 

SOC 

0% 

20% 

40% 

60% 

80% 

100% 

—1  °C 

o°c 

25  °C 

45  °C 

K(U) 

K  2.1) 

R(  1.2) 

K(4,6) 

(7) 


Side-reaction  current,  Js  adversely  affects  cell  performance  by 
generating  mixed  potential  in  the  present  model.  It  is  added  to  the 
model  in  order  to  consider  cell  life  degradation.  Although  an  actual 
lithium-ion  cell  body  is  a  composite  structure  made  of  several 
layers  such  as  electrodes,  separator  and  current  collector,  it  is 
assumed  that  the  dominant  electrical  current  flow  occurs  in  the  in¬ 
plane  direction  along  current  collector  which  has  the  lowest  elec¬ 
tric  resistance  in  the  present  model.  Therefore,  the  average  electric 
conductivity  of  the  cell  body  domain  in  Eq.  (6)  is  acquired  by  scale 
factors  as  follows: 


Vp 


ap  •  oM  where  ap  =  nAi-f°jyHM-foii 


an  =  an-  aCu  where  an 


^Cu-foil '  ffcu-foil 

Ffcell 


(8) 


This  is  because  the  present  multi-D  model  is  based  on  ECM,  and 
it  is  not  necessary  to  consider  the  detailed  structure  of  each  cell 
layer.  In  the  positive  tab  which  is  made  of  aluminum,  the  following 
electric  equation  is  solved. 


(TMV2<Pp  =  0 

In  the  same  manner,  the  negative  tab  equation  is 


(9) 
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Fig.  5.  Schematic  of  muiti-D  lithium-ion  cell  model  which  has  constituting  ECMs  in  each  mesh  node:  (left)  uni-directional  cell,  (right)  bi-directional  cell. 


Table  2 

Specification  of  the  lithium-ion  cell  model  for  the  present  study. 


Parameter 

Sym. 

Value 

Unit 

Cell  length 

f-cell 

190 

mm 

Cell  width 

Wcell 

142 

mm 

Cell  thickness 

Ffcell 

5.5 

mm 

Thickness  of  negative  current  collector 

HCu 

20 

pm 

Thickness  of  positive  current  collector 

Hai 

20 

pm 

Thickness  of  negative  electrode 

% 

50 

pm 

Thickness  of  positive  electrode 

55 

pm 

Thickness  of  separator 

^sep 

20 

pm 

Lead-tab  length 

f-tab 

40 

mm 

Lead-tab  width 

Wtab 

40 

mm 

Lead-tab  thickness 

Htab 

0.2 

mm 

Nominal  capacity  (1C  capacity) 

cap0 

15.0 

Ah 

Cutoff  voltage  (max) 

VW 

4.15 

V 

Cutoff  voltage  (min) 

^min 

3.0 

V 

Density  of  cell 

P 

2400 

kg  ITT3 

Specific  heat  of  cell 

cp 

1015 

J  kg'1  K1 

Fig.  6.  Experimental  validation  of  the  multi-D  model  according  to  continuous 
discharging. 


Fig.  7.  Comparison  of  cell  voltage  profile  between  simulation  and  experiment  under 
repetitive  80  A  charging/discharging  current  pulse  (SOQnit  =  60%) :  (a)  Isothermal 
model,  (b)  non-isothermal  model. 
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Fig.  8.  Comparison  of  ceil  temperature  (point  B)  between  simulation  and  experiment 
under  repetitive  80  A  charging/discharging  current  pulse  (Tamb  =  288  K, 
h  =  0.18  W  m"2  K1). 


Table  3 

Experimental  validation  summary  (US06  driving  test). 


T(°C) 

Discharging  &  rest 

US06  (5  cycles) 

Avg.  error 

(%) 

Std.  deviation 

(%) 

Avg.  error 

(%) 

Std.  deviation 

(%) 

45 

0.4 

0.4 

0.4 

0.4 

25 

0.5 

0.4 

0.7 

0.5 

15 

- 

- 

1.5 

1.1 

5 

- 

- 

2.4 

1.9 

-10 

1.2 

0.8 

2.9 

1.8 

<7cuV2<2>„  =  0  (10) 

Note  the  above  electrical  equations  came  from  the  Maxwell 
equation  by  ignoring  magnetic  field  effects.  In  the  cell  body  domain 
where  electrochemical  reaction  occurs  (J  ¥=  0),  local  cell  voltage, 
Vceii (x,y,z)  from  ECM  Eq.  (1)  is  same  as  the  potential  difference 
calculated  from  Eq.  (6)  as  follows: 

Vceii (x,y,z)  =  U(x,y,z)  -  Vs(x,y,z)  -  V1(x,y,z)  -  V2(x,y,z) 

=  <Pp(x,y,z)  -  <Pn(x,y,z) 

(11) 

From  Eq.  (11),  Vs  can  be  defined  as  follows: 

Vs  =  U-{0p-<Pn}-V1-V2  (12) 


Fig.  9.  Experimental  validation  of  the  present  multi-D  model:  (a)  US06  driving  test  @25  °C,  (b)  US06  driving  test  @15  °C,  (c)  US06  driving  test  @-10  °C. 
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0.3044E+06 
0.3039E+06 
0.3034E+06 
0.3029E+06 
0.3024E+06 
0.301 9E+06 
0.301 4E+06 
0.3009E+06 
0.3004E+06 
0.2999E+06 
0.2995E+06 
0.2990E+06 
0.2985E+06 
0.2980E+06 
0.2975E+06 


0.5240 

0.5230 

0.5220 

0.5210 

0.5200 

0.5190 

0.5180 

0.5170 

0.5160 

0.5150 

0.5140 

0.5130 

0.5120 

0.5110 

0.5100 


(c) 


3.789 

3.788 

3.788 

3.787 

3.786 

3.785 

3.784 

3.784 

3.783 

3.782 

3.781 

3.780 

3.780 

3.779 

3.778 


Fig.  10.  Distribution  contour  of  (a)  volumetric  current  density  distribution  (unit:  A  m  3),  (b)  SOC,  and  (c)  OCV  under  3C  rate  discharging  condition  @600  s. 


Therefore,  volumetric  current  density  at  point  (x,y,z)  can  be 
found  from  Eq.  (2)  as 


In  the  above  energy  equation,  thermal  conductivity  of  the  cell 
body  is  acquired  from 


Jrxn 


VceirSs 


(13) 


As  the  ECM  parameters  are  strongly  affected  by  temperature,  it 
is  necessary  to  couple  the  energy  equation  as  follows: 


pep  f  =  V(kVT)  —  h(T  —  Too)  +  q" 

where  q'"  =  q’’’xn  +  <j”hm  =  j(u  -  V  +  T^fj  +  JM^)2 

(14) 


kin-plane  =  £  =  30(W  m-1  I<  ') 

i 

l  x-1  (15) 

kthru-plane  =  fr^E&j  =  0.8(W  ITT1  K_1) 

It  is  noteworthy  that  heat  conduction  occurs  faster  in  the  in¬ 
plane  direction  than  in  the  thru-plane  direction.  The  reaction 
term  q"  in  Eq.  (14)  becomes  zero  in  the  positive  tab  and  negative 
tab  domain. 
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4.878 

4.877 

4.876 

4.875 

4.873 

4.872 

4.871 

4.870 

4.869 

4.868 

4.867 

4.866 

4.865 

4.864 

4.862 


0.9998 

0.9987 

0.9976 

0.9965 

0.9954 

0.9942 

0.9931 

0.9920 

0.9909 

0.9898 

0.9887 

0.9876 

0.9865 

0.9854 

0.9842 


Fig.  11.  Contour  of  positive  and  negative  potential  distribution  (3C  discharging)  at  100  s  (unit:  V):  (a)  Positive  potential,  (b)  Negative  potential. 


Boundary  conditions  for  galvanostatic  operation  is  given  as 
follows: 


-a 


eff^s 

Cu  ax 


Cu-tab 


eff  d(Ps 


UA\ 


OX 


Al-tab 


I 

^tab 


(16) 


The  present  model  can  also  simulate  other  cell  operation  modes 
such  as  potentiostatic,  power-driven,  and  pulse  mode  by  changing 
boundary  conditions  above. 


2.3.  Degradation  model 

Battery  degradation  is  expressed  as  capacity  fade  and  imped¬ 
ance  growth.  Degradation  can  be  caused  when  charging  and 


Table  4 

Degradation  simulation  cases. 


Case 

no. 

Run  mode 

C-rate 

T 

(°C) 

k 

(W  ITT1  K1 

Tab  position 

) 

Remark 

1 

100  cycles 

3 

25 

30.0 

Uni-directional 

Baseline  case 

2 

100  cycles 

5 

25 

30.0 

Uni-directional 

Effect  of  C-rate 

3 

100  cycles 

3 

45 

30.0 

Uni-directional 

Effect  of 
temperature 

4 

100  cycles 

3 

25 

5.0 

Uni-directional 

Effect  of 
thermal 
conductivity 

5 

100  cycles 

3 

25 

30.0 

Bi-directional 

Effect  of 
lead-tab 
position 

307.2 
307.0 

306.8 

306.5 

306.3 

306.1 

305.9 

305.6 

305.4 

305.2 

304.9 

304.7 


304.5 

304.3 


304.0 


320.6 

319.7 

318.7 

317.8 

316.8 

315.9 
315.0 
314.0 

313.1 

312.1 

311.2 

310.2 

309.3 

308.3 

307.4 


Fig.  12.  Contour  of  cell  temperature  distribution  (unit:  K):  (a)  3C  rate  @600  s,  (b)  5C  rate  @300  s. 
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discharging  of  a  battery  is  repetitively  cycled,  which  is  called 
cycle-life  loss.  When  a  battery  is  stored  in  a  place  for  long  time 
without  operation,  degradation  may  also  occur,  which  is  called 
calendar-life  loss.  Usually,  both  degradation  modes  take  place  in  a 
combined  manner  during  the  entire  life  time  of  battery.  In  order  to 
develop  long-lasting  batteries,  it  is  important  to  understand  the 
cause  and  mechanism  of  battery  degradation.  Degradation  origi¬ 
nates  from  several  causes  as  lithium-ion  battery  is  a  very 
complicated  system  [12].  It  is  still  unclear  how  causes  of  degra¬ 
dation  are  related  and  how  they  affect  the  cell  life.  Therefore,  it  is 
not  easy  to  construct  a  physics-based  degradation  model  although 
several  studies  have  been  conducted  under  restricted  conditions 
[13-17].  Empirical  models  [10  can  be  constructed  relatively  easily 
although  it  is  not  easy  to  collect  experimental  data  for  long  period. 
Semi-empirical  models  [11]  can  be  applied  to  more  general 
simulation  cases  than  empirical  models  because  model  parame¬ 
ters  are  flexible. 

Recently,  the  sizes  of  lithium-ion  cells  have  become  larger.  In 
this  case,  local  degradation  distribution  becomes  more  important, 
which  may  comes  from  cell  size,  shape  or  lead-tab  position.  As  most 
of  degradation  model  listed  above  are  1-D  or  lumped-scale  model, 
they  cannot  predict  local  degradation  coming  from  dimensional 
effect.  Multi-D  cell  model  has  capability  of  predicting  local  degra¬ 
dation  distribution.  Various  empirical  or  semi-empirical  degrada¬ 
tion  models  can  be  incorporated  with  the  present  multi-D  model.  In 
this  paper,  a  simple  empirical  degradation  model  which  predicts 
capacity  fade  and  resistance  growth  according  to  cell  operation  and 
exposure  to  temperature  is  presented  to  show  local  degradation 
distribution. 

The  present  model  assumes  nominal  capacity  at  location  (x,  y,  z) 
is  a  function  of  total  accumulated  coulomb  (Q)  as  proposed  by 
Bloom  et  al.  [10]. 

cap  =  cap0(^i  -  a-Q}/2^  (17) 

Eq.  (17)  imitates  cell  capacity  fade  by  available  lithium  loss.  Local 
SOC  is  determined  based  on  the  updated  nominal  capacity  in 
Eq.  2(17). 


soc  =  soc0  -  =L_-  J  J(t)  dt  (18) 

time 

Accumulated  coulomb  is  sum  of  reaction  current  C/rxn)  and  side- 
reaction  current  (Js ).  Side  reaction  current  is  added  to  consider 
calendar  life. 


Q«)  -  /  (/rxn  +Js)  dt 

t 

State  of  health  is  defined  as 


(19) 


soh 


cap 
^ Po 


(20) 


Three  resistance  terms  are  modified  as  a  function  of  SOH  to 
consider  resistance  growth. 


Rs  =  Rs,0{l+Ml-soh)} 

R,  =Rh  0{l+/V(l-soh)}  (21) 

R2  =R2, 0^+ Ml -  soh)} 

In  Eq.  (21 ),  /3  is  adjustable  fitting  parameter.  Because  the  present 
model  solves  governing  equation  at  each  time  step  with  given 
boundary  condition,  the  model  can  simulate  cell  degradation  under 
various  operating  conditions  under  any  current  (or  power)  load 
profile  and  thermal  profile.  Although  relatively  simple  degradation 
model  for  multi-D  application  is  introduced  in  this  study,  more 


sophisticated  life  models  18]  can  be  implemented  into  the  present 
multi-D  model. 

3.  Result  and  analysis 

Specification  of  a  15  Ah-class  Li-ion  cell  used  for  the  present 
study  is  presented  in  Table  2.  The  lithium-ion  cell  has  a  positive 
electrode  made  of  layered  NMC  (Lithium  nickel  manganese  cobalt) 
and  manganese  spinel  while  graphite  and  carbon  are  blended  to 
construct  a  negative  electrode.  The  multi-D  model  is  solved  by 
finite  volume  method  (FVM).  Solution  is  assumed  to  converge 
when  residuals  reaches  1CT9. 

3.1  Baseline  multi- dimensional  cell  model 

Experimental  validation  of  the  present  multi-D  model  under 
several  C-rate  discharging  conditions  is  presented  in  Fig.  6.  The 
model  shows  pretty  good  accuracy  up  to  7C  discharging  condition. 
However,  simulation  result  starts  to  deviate  from  the  experimental 
data  above  7C-rate.  This  is  because  it  is  difficult  to  handle  mass 
transport  limiting  phenomena  with  standard  ECM.  Additional  var¬ 
iables  such  as  Warburg  circuit  will  be  helpful  to  resolve  this  prob¬ 
lem,  which  is  not  covered  in  this  study. 

In  order  to  validate  dynamic  voltage  response  of  the  model 
against  experimental  data,  a  repetitive  charging/discharging  cur¬ 
rent  pulse  test  during  2000  s  is  conducted  as  shown  in  Fig.  7(a). 
When  temperature  is  maintained  constant,  simulation  result  shows 
almost  same  magnitude  of  voltage  response  since  same  amount  of 
coulomb  is  charged  and  discharged  during  the  test.  However,  the 
actual  cell  shows  slightly  diminishing  magnitude  of  voltage 
response.  This  is  because  the  cell  resistance  decreases  as  the  cell 
temperature  increases  during  the  pulse  test.  Temperature  profile 
during  the  pulse  test  is  shown  in  Fig.  8.  When  the  thermal  model  is 
applied  to  the  test  scenario,  ECM  parameters  are  updated  according 
to  the  temperature  and  the  voltage  response  of  the  model  matches 
well  with  the  experimental  data  as  shown  in  Fig.  7  (b).  Driving 
simulation  according  to  US06  driving  profile  is  conducted  at  tem¬ 
perature  of  25  °C,  15  °C,  and  -10  °C  (see  Fig.  9).  Note  15  °C 
and  -10  °C  were  not  included  in  Table  1,  which  means  these  two 
simulation  cases  are  conducted  with  bi-linearly  interpolated  pa¬ 
rameters.  Although  interpolated  model  parameters  are  used, 


Fig.  13.  Simulated  voltage  profile  of  the  cell  for  the  first  cycle  (case  #1). 
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Fig.  14.  SOH  according  to  cyclic  operation  mode  of  the  multi-D  model:  (a)  case  #1,  (b)  case  #2,  (c)  case  #3,  (d)  case  #4,  (e)  case  #5. 


simulation  result  agrees  well  with  experiment  data.  Experimental 
validation  under  other  conditions  is  summarized  in  Table  3.  Model 
accuracy  slightly  decreases  when  the  cell  operates  at  low  temper¬ 
ature.  This  is  because  model  parameters  become  more  sensitive  to 
the  ambient  temperature  under  low  temperature  conditions.  In 
addition,  it  is  difficult  to  get  accurate  model  parameters  at  low 
temperature  because  maintaining  isothermal  condition  of  the 
operating  cell  is  not  easy. 

Volumetric  current  density  (reaction  rate)  during  3C  rate  dis¬ 
charging  condition  is  presented  in  Fig.  10(a).  It  is  observed  that 
current  density  is  strongest  near  lead-tab  because  electric  current 
produced  in  the  electrode  is  collected  by  lead-tab.  As  the  current 
density  is  high  near  lead-tab  region,  local  SOC  is  low  there  also  as 
shown  in  Fig.  10(b).  Local  OCV  can  be  found  according  to  local  SOC 


(Fig.  10(c).  OCV  near  lead-tab  drops  faster  than  other  region  since 
reaction  rate  is  high  there  during  discharging.  These  results  agree 
with  our  previous  intuitive  expectations.  Potential  field  in  each 
electrode  is  shown  in  Fig.  11.  Potential  field  distribution  is  strongly 
affected  by  the  effective  electric  conductivity  in  Eq.  (8).  When  the 
current  collector  is  thick  and  made  of  more  conductive  material,  the 
effective  electric  conductivity  becomes  high.  As  a  result,  it  improves 
the  uniformity  of  potential,  OCV,  and  SOC  distribution  in  each 
electrode,  which  results  in  reduced  voltage  loss.  Amount  of  local 
heat  generation  in  the  cell  is  calculated  from  Eq.  (14)  with  the 
current  density  and  potential  distribution  that  are  calculated.  Then, 
temperature  distribution  is  obtained  by  solving  energy  equation.  As 
the  present  cell  model  is  thermally  coupled,  local  temperature  is 
used  to  find  appropriate  temperature-dependent  model 


S.  Jung,  D.  Kang  /  Journal  of  Power  Sources  248  (2014)  498-509 


507 


0.9460 

0.9459 

0.9458 

0.9457 

0.9456 

0.9455 

0.9454 

0.9453 

09452 

0.9450 

0.9449 

0.9448 

0.9447 

0.9446 

0.9445 


0.9460 

0.9459 

0.9458 

0.9457 

0.9456 

09455 

09454 

09453 

09452 

0.9450 

0.9449 

0.9448 

09447 

0.9446 

0.9445 


0.9429 

0.9429 

0.9428 

0.9427 

0.9427 

0.9426 

0.9425 

0.9425 

0.9424 

0.9423 

0.9423 

0.9422 

0.9421 

0.9421 

0.9420 


Fig.  15.  Contours  of  SOH  during  cyclic  degradation  test  (after  100  cycles):  (a)  case  #1,  (b)  case  #4,  (c)  case  #5. 


parameters  such  as  resistances.  Fig.  12  shows  temperature  distri¬ 
bution  of  the  cell.  As  the  electric  current  converges  to  lead-tab, 
temperature  is  highest  near  lead-tab  region.  Temperature  distri¬ 
bution  in  the  cell  body  domain  is  strongly  affected  by  thermal 
conductivity.  If  thermal  conductivity  is  too  low,  temperature  dif¬ 
ference  in  a  cell  becomes  large  and  it  may  deteriorate  state  of  heath 
in  a  cell. 

3.2.  Degradation  effect 

In  this  part,  the  effect  of  degradation  on  SOH  of  a  cell  is  dis¬ 
cussed  with  the  degradation  model.  Table  4  lists  simulation  cases 
conducted  to  examine  degradation  effect  on  the  cell  performance 
according  to  C-rate,  temperature,  thermal  conductivity,  and  lead- 
tab  position. 

In  order  to  conduct  a  cyclic  degradation  simulation,  a  single 
cycle  with  the  baseline  cell  model  is  prepared  as  shown  in  Fig.  13. 
First,  the  cell  discharges  until  the  cell  voltage  reaches  lower  cutoff 
voltage.  After  taking  rest  for  0.5  h,  the  cell  is  charged  until  the  cell 
voltage  reaches  upper  cutoff  voltage.  Finally,  the  cell  is  charged 
with  constant  voltage  until  charging  current  reaches  5%  of  1C 
current.  This  run  is  cycled  100  times  to  examine  the  degradation  of 
the  cell. 


Average  SOH  of  the  baseline  cell  according  to  cycle  number  is 
presented  in  Fig.  14(a)— (c).  It  is  observed  that  SOH  decreases  faster 
when  the  cell  operates  under  severe  operating  condition  such  as 
high  C-rate  and  high  temperature.  However,  calculation  result 
shows  that  standard  deviation  (sigma)  of  SOH  is  low  under  severe 
operating  condition.  In  other  words,  uniformity  of  SOH  becomes 
improved.  When  the  cell  operates  under  severe  condition,  cell 
temperature  is  elevated  and  the  cell  resistances  are  reduced  ac¬ 
cording  to  high  temperature.  As  the  cell  resistances  are  reduced, 
uniformity  of  local  overpotential  in  the  cell  is  improved  resulting  in 
uniform  accumulated  coulomb  i.e.  state  of  health.  This  result  hap¬ 
pens  since  the  present  model  parameters  are  strongly  affected  by 
temperature.  If  thermal  conductivity  of  a  cell  is  low,  temperature 
uniformity  will  be  worsened  which  leads  to  low  SOH  uniformity  as 
shown  in  Fig.  14(d).  This  result  implies  high  thermal  conductivity  of 
a  cell  mitigates  local  degradation  due  to  improved  uniformity  of 
temperature  distribution.  Calculation  result  revealed  that  the  cell 
with  bi-directional  lead-tabs  shows  the  best  degradation  perfor¬ 
mance  (the  highest  remaining  SOH  and  the  highest  SOH  unifor¬ 
mity)  compared  to  other  cases.  Local  distribution  contour  of  SOH  is 
shown  in  Fig.  15.  SOH  in  the  region  near  lead-tab  is  lower  than  other 
region  since  the  region  near  lead-tab  experienced  relatively  severe 
operation  (high  current  density  and  high  temperature)  compared 
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Fig.  16.  Cell  resistance  (R2)  variation  according  to  cyclic  operation  case  #1  (small  graph 
shows  cell  resistance  variation  in  a  single  cycle:  discharging/rest/charging/rest. 


to  the  other  regions.  This  local  degradation  will  be  mitigated  if  lead- 
tab  is  wide  and  thick,  which  may  lower  heat  generation  and  current 
density  near  lead-tab.  The  cell  with  bi-directional  lead-tabs  shows 
better  distribution  of  temperature  and  SOH  than  the  cell  with  uni¬ 
directional  lead-tabs  does.  In  the  cell  with  bi-directional  lead-tabs, 
temperature  near  positive  lead-tab  had  been  maintained  relatively 
higher  than  other  regions,  which  results  in  low  SOH  there  (see 
Fig.  15(c)). 

Cell  resistances  increase  as  the  cell  degradation  goes  on.  Fig.  16 
shows  variation  of  a  cell-averaged  ECM  resistance  when  the  cell 
operates  for  100  cycles.  Magnified  figure  in  the  plot  shows  resis¬ 
tance  variation  in  a  single  cycle.  Note  resistance  difference  between 
SOC  0%  and  SOC  100%  in  a  cycle  increases  as  the  cell  degradation 
goes.  As  the  cell  resistance  increases,  overpotential  (U—V)  increases 
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Fig.  17.  Temperature  at  the  center  of  the  cell  during  cyclic  operation  case  1  (small 
graph  shows  temperature  variation  in  a  single  cycle:  discharging/rest/charging/rest). 


and  the  amount  of  heat  generation  increases  also  according  to  Eq. 
(14).  This  leads  to  temperature  increase  as  SOH  of  the  cell  is  worsen 
as  shown  in  Fig.  17.  Therefore,  it  is  important  to  consider  temper¬ 
ature  increase  due  to  cell  degradation  when  designing  battery  pack 
cooling  system. 

4.  Conclusion 

A  multi-D  model  of  lithium-ion  cell  has  been  developed.  The 
model  is  based  on  ECM  and  it  is  capable  of  simulating  various  cell 
operation  modes  with  good  accuracy.  In  order  to  consider  effect  of 
temperature  and  SOC  on  the  cell  performance,  model  parameters 
were  bilinearly  interpolated.  Simulation  results  show  good  agree¬ 
ment  with  experiment  data  under  various  operation  conditions.  Local 
degradation  effects  have  been  elucidated  by  applying  an  empirical 
degradation  model  to  the  present  multi-D  model.  Degradation 
became  worsen  under  severe  condition  such  as  high  storage  tem¬ 
perature  and  high  C-rate.  It  was  observed  that  low  thermal  conduc¬ 
tivity  of  a  cell  is  harmful  to  uniformity  of  cell  SOH.  Local  distributions 
of  potential,  SOC,  SOH  have  been  examined,  which  is  difficult  to 
observe  with  conventional  experimental  methods. 

The  present  model  has  potential  to  be  extended  to  various  ap¬ 
plications.  For  example,  by  separating  positive  and  negative  elec¬ 
trode,  the  present  multi-D  model  can  be  extended  to  study  battery 
safety  such  as  external  shorts  or  nail-penetration.  In  addition,  the 
developed  model  can  be  extended  to  transient  CFD  analysis  of 
large-scale  battery  pack,  which  is  necessary  for  designing  battery 
pack  cooling  system. 

Through  an  extensive  work  of  model  validation,  it  was 
confirmed  that  the  present  multi-D  model  works  well  with  single¬ 
phase  active  materials  such  as  NMC  (Lithium  nickel  manganese 
cobalt),  LCO  (Lithium  cobalt  oxide),  LMO  (Lithium  manganese 
spinel),  and  graphite.  But,  ECM  sub-model  inserted  into  the  present 
multi-D  model  needs  to  be  modified  if  the  electrode  is  blended 
with  phase-changing  materials  such  as  LiFeP04,  and  LTO  (Lithium 
titanate  oxide)  in  order  to  reflect  phase-transition  phenomena. 

List  of  symbols 


A  area  (cm2) 

cp  specific  heat  (J  kg-1  I<-1) 

C  capacitance  (Farad) 

cap  capacity  (Ah) 

Ea  activation  energy  (J  mol-1) 

F  Faraday  constant  (96487  C  mol-1) 

h  heat  transfer  coefficient  (W  m-2  K-1) 

H  thickness  (cm) 

I  current  (A) 

J  volumetric  current  density  (A  cm-3) 

k  thermal  conductivity  (W  m-1  K-1) 

L  length  (cm) 

m  mass  (kg) 

n  number 

R  resistance  (ohm),  universal  gas  constant 

(8.314  J  mol-1  K-1) 
soc  state  of  charge 

soh  state  of  health 

t  time  (sec) 

T  temperature  (I<) 

V  volume  (cm3) 

w  width  (cm) 

U  open  circuit  potential  (V) 

q  heat  source  (W) 

Q.  charge  (Coulomb) 
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V  voltage  (V) 

V  volume  (m3) 

Greek 

a  data  fitting  coefficient 

7]  overpotential  (V) 

p  density  (kg  m-3) 

a  electric  conductivity  (S  cm-1) 

0  potential  (V) 


Superscripts  and  subscripts 

0 

initial 

amb 

ambient 

cell 

cell 

eff 

effective 

n 

negative 

P 

positive 

ref 

reference 

rxn 

reaction 

s 

serial,  side  reaction 

sep 

separator 

tab 

lead-tab 
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